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Department of Mechanical and Aerospace Engineering 
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Abstract 

The purpose of this research is to continue our efforts in advancing the state of 
knowledge in large eddy simulation (LES), direct numerical simulation (DNS) and 
Reynolds averaged Navier Stokes (RANS) methods for the computational analysis of 
high-speed reacting turbulent flows. In the second phase of this work, covering the 
period: September 1, 1993 - September 1, 1994, we have focused our efforts on two 
research problems: (1) Developments of “algebraic” moment closures for statistical 
descriptions of nonpremixed reacting systems, (2) Assessments of the Dirichlet fre- 
quency in presumed scalar probability density function (PDF) methods in stochastic 
description of turbulent reacting flows. This report provides a complete description of 
our efforts during this past year as supported by the NASA Langley Research Center 
under Grant NAG-1-1122. 


Technical Monitor: 


Dr. J. Philip Drummond, Hypersonic Propulsion Branch, Tel: 804-864-2298 is the Technical 
Monitor of this Grant. 
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1 Introduction 


While this program is associated with research on several problems of current interest in 
reacting turbulent flows, its primary objective is to facilitate the use of LES for the com- 
putational analyses of practical high speed transport. With the evolution of our work, 
sponsored by NASA for the past four years, it has become clear that to meet the needs of 
this program it is required to make use of advanced, yet practical, stochastic and statistical 
procedures in conjunctions with reliable computational procedures. The importance of the 
these methods have become clear through the results of our most recent work on LES of 
turbulent reacting flows . 1 It is now well-established that many of the conventional closures 
which work reasonably well in RANS 2 fail in LES. Thus, it has become clear that the first 
step to be taken in any meaningful LES of reacting flows is to make sure that the moments 
“up to the second order level” are modeled accurately. It must be noted that in none of the 
previous contributions in LES, was this the subject of a detailed study (see Ref. 3 ). In fact, 
in almost all previous and concurrent contributions, only the “first” subgrid scale moment of 
the transport variables have been the subject of modeling. That is, the approach is based on 
a variation of the Smagorinsky 4 closure. Unfortunately the application of such closures does 
not work in LES of reacting flows. In our efforts within the past year, we have made use of 
the recent work of one of the Co-PI’s of this proposal (D.B. Taulbee) in developing improved 
algebraic models. In his previous work, Taulbee 5 shows that with the modeled dynamics 
equations for the Reynolds stresses, it is possible to develop an improved explicit algebraic 
Reynolds stress model that can predict many of the flow features more accurately than the 
conventional models. This improvement is due to the fundamentals of the approach in that 
the closure is based on transport equations for the higher-order moments. Therefore, more 
physics is embedded in the equations. Furthermore, the extra degree of freedom provided by 
the closure allows more adaptability in its optimization for predictive analysis. In fact, some 
sample results in Refs . 1,6 do show that these new models are superior to previous closures 
in LES, even though the difference is not significant in RANS (also see Ref. 2 ). 

In LES of reacting flow phenomena in addition to the first two moments, the PDF of the 
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scalars within the subgrid must also be specified. Thus, we have continued our work on 
PDF modeling as well. Based on our extensive work in the first phase of this research 7- ® 
we have established that the most meaningful means of utilizing the PDF in LES is by 
mean of describing and solving the transport equations governing its evolution. However, 
it has also become clear to us the most recent PDF procedure based on the Amplitude 
Mapping Closure of Kraichnan 10 - 11 cannot be used for any type of practical applications 
(nor can they be utilized in basic applications in flows with nonequilibrium chemistry). 
Therefore, we have decided to rely and focus on the conventional Coalescence/Dispersion 
methodology. 12-14 Work is underway in developing numerical schemes based on the Monte 
Carlo methods to solve the PDF within the subgrid in LES of reacting flows. We do not have 
any substantiated results to report at this time. However, we do have extensive new results 
in PDF modeling based on “presumed methods”. In particular, we have made a rather 
extensive study in determining the capabilities and drawbacks of the Dirichlet frequency in 
probability modeling of nonpremixed reacting flows. Our reason for this study is due to the 
interest at NASA LaRC in the approach. In particular, we cite the work of Gaffney et a/. 15 ' 16 
who have used the Dirichlet density for modeling of high speed reacting systems. We feel 
that a detailed appraisal of this versatile density is in order, and we have made use of all 
recent laboratory data for this purpose. 


2 Summary of Achievements 

Appendix I and Appendix II provide a complete write-up on our accomplishments on the two 
components of this program. For the convenience of the reader, here a summary is provided: 

Our efforts in (1) have been devoted towards developing closures which can be used for mod- 
eling of the “second order moments” in the contexts of both RANS and LES. In particular, 
explicit algebraic scalar flux models that are valid for three-dimensional turbulent flows are 
derived from a hierarchy of second-order moment closures. The mathematical procedure 
is based on the Cayley- Hamilton theorem and is an extension of the scheme developed by 
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Pope 17 and Taulbee. 5 Several closures for the pressure-scalar gradient correlations are con- 
sidered and explicit algebraic relations are provided of the velocity-scalar correlations in 
both non-reacting and reacting flows. In the latter, the role of the Damkohler number is ex- 
plicitly exhibited in turbulent flows with nonpremixed reactants. The relationship between 
these closures and traditional models based on the linear gradient diffusion approximation is 
theoretically established. The results of model predictions are assessed via comparison with 
available laboratory data. 

In efforts related to (2), the recent experimental data in Ref. 18 pertaining to compositional 
structure of a turbulent reactive scalar mixing layer are reproduced by a mathematical- 
computational procedure utilizing the Pearson family (PF) of univariate and multivariate 
PDF’s. By detailed comparisons against these data and some additional data generated by 
DNS of a spatially developing reacting mixing layer, an appraisal is made of the applicability 
and the extent of validity of PF for statistical description of reactant fluctuations. In accord 
with the experiment, a chemical reaction of the type A + B —* Products is simulated in 
isothermal, incompressible flows. A wide range of the Damkohler number is considered 
including both frozen and equilibrium chemistry flows. The comparison of the results with 
laboratory data indicates that PF generated PDF’s are very convenient, in the absence of 
better alternatives, for modeling the influence of turbulence on the reactant conversion rate. 
In particular, the Dirichlet frequency provides the most reasonable means of portraying 
the multivariate scalar PDF. The extent of agreement improves as the magnitude of the 
Damkohler number is increased. A more detailed comparative assessment of the model 
predictions against DNS data confirms the relative superiority of the Dirichlet PDF even 
though it is mathematically shown that this frequency is invalid in equilibrium flows. With 
the use of the PF generated PDF, the autospectral density function and the cross-spectra 
density function (including both the coherence and the phase) of the reacting scalars under 
equilibrium chemistry are related to frequency spectra of a conserved mixture fraction. This 
relation is very convenient and is favored over previous models for predicting the spectral 
characteristics of reacting scalars in the central region of the laboratory mixing layer and in 
any other homogeneous flow configurations. 
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The results of our work in both (1) and (2) are complete. Two papers have been written 
and are submitted for publication. 


3 Work in Progress 

In addition to the main two topics discussed in the appendixes, work is in progress in two 
other research areas: (1) LES by means of the Monte Carlo solution of the subgrid scalar 
PDF, and (2) LES by means of mechanistic models. In (1), we have just completed the 
development of a combined finite difference-Monte Carlo procedure for solving the PDF of 
scalar variables within the subgrid. At this point, we are faced with one major problem and 
that is the Monte Carlo solution of the advection term in the PDF transport. We have used 
a first order upwinding scheme for this procedure. This methods is stable, but introduces a 
significant amount of artificial diffusivity. In fact the amount of diffusivity introduced by the 
scheme is more that by the subgrid closure. Unfortunately, none of the higher order schemes, 
including all of those in Refs. 19,20 are amenable to Monte Carlo discretization. This serious 
problem is the subject of our current investigation. 

In (2), due to our recent progress in the area of diffusion flamelet modeling (DFM) 21 23 of 
nonpremixed reacting flows, we have devoted a small portion of our time toward exploring 
the possibility of developing subgrid scale closures by means of DFM. Based on our earlier 
findings reported in Ref., 24 we expect that the model should work reasonably well in the 
flamelet region; that is when the chemistry is fast and the thickness of the flame is very small. 
In such cases, the LES methodology based on DFM is expected to be satisfactory. For this 
purpose we have initiated a task in which typical mixing controlled chemistry models such 
as those developed in Ref. 25 are used in a dynamic subgrid model 26 for LES of homogeneous 
reacting flows. We have completed the mathematical formulation to accomplish this task; 
but we do not have any significant numerical data to report at the present time. 

We hope to have results in each of these two endeavors before the deadline for our next 
semi-annual report. 
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4 Personnel 


Dr. Peyman Givi is as the PI of this project. Drs. and Dale B. Taulbee and Cyrus K. Madnia 
serve as the Co-PI’s. Of course, as before, Dr. Givi is responsible for a timely and successful 
completion of all the tasks. 

There are two Ph.D. candidates who are supported by this grant: Mr. Virgil Adumitroaie and 
Mr. Craig Steinberger. In addition, the following students have contributed to various parts 
of the project, but are not financially supported by the NASA grant: Mr. Geoffrey Shieh, 
Mr. Sean Garrick and Mr. George Sabini. All these students have recently completed the 
requirements for M.S. degree at SUNY-Buffalo. Mr. Shieh and Mr. Garrick are continuing 
their studies toward Ph.D. at SUNY-Buffalo and are expected to continue the the work of 
Virgil Adumitroaie and Craig Steinberger. The latter two are expected to complete the 
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APPENDIX I 


Explicit Algebraic Scalar-Flux Models for 
Turbulent Reacting Flows 


V. Adumitroaie, D.B. Taulbee and P. Givi 
Department of Mechanical and Aerospace Engineering 
State University of New York at Buffalo 
Buffalo, NY 14260-4400 


Abstract 

Explicit algebraic scalar flux models that are valid for three-dimensional turbulent 
flows are derived from a hierarchy of second-order moment closures. The mathemat- 
ical procedure is based on the Cayley- Hamilton theorem and is an extension of the 
scheme developed by Pope 1 and Taulbee. 2 Several closures for the pressure-scalar gra- 
dient correlations are considered and explicit algebraic relations are provided of the 
velocity-scalar correlations in both non-reacting and reacting flows. In the latter, the 
role of the Damkohler number is explicitly exhibited in turbulent flows with non- 
premixed reactants. The relationship between these closures and traditional models 
based on the linear gradient diffusion approximation is theoretically established. The 
results of model predictions are assessed via comparison with available laboratory data. 

PACS: 47.27. Eq (Turbulence simulation and modeling), 47.70.Fw (Chemically re- 
active flows), 47.27. Qb (Turbulent diffusion), 47.27. Nz (Boundary layer and shear tur- 
bulence), 47.27. Wg (Jets) 02.10.Sp (Matrix theory). 


I Introduction 

Despite extensive recent contributions in direct and large eddy simulations of turbulent 
reacting flows, the application of such simulations is limited to “simple flows”. 3-5 Based on 
this fact, it is now widely recognized that the “statistical” approach is still the most practical 
means in computational turbulence, and future capabilities in predictions of engineering 
turbulent combustion systems depend on the extent of developments in statistical modeling. 

The literature on computational prediction of nonreactive turbulent transport- is abun- 
dant with schemes based on single-point statistical closures for moments up to the second- 
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order. 6-10 The presence of scalar contaminant and/or chemical reactions generates additional 
length and time scales which have to be considered. 3,11-16 To account for these scales in a 
second-order moment formulation the solution of a large number of transport equations is 
required. This could potentially make the approach less attractive, but can be alleviated 
by utilizing “algebraic closures”. 1,2,17-22 Such closures are either derived directly from the 
modeled transport equations of the respective moments, or other types of representations 
that lead to anisotropic “eddy-diffusivities”. 23-26 In this manner, the number of equations is 
reduced but the accuracy of the second-order formulation and the versatility of the approach 
is preserved. 

In this work, we expand upon the formulation developed by Taulbee 2 (also see Ref. 22 ) to pro- 
vide explicit algebraic relations for the flux of scalar variables in isothermal turbulent flows. 
Both nonreacting and reacting flows are considered. In the latter, a second-order, irreversible 
chemical reaction of the type A + B — » P is considered in turbulent flows with initially segre- 
gated reactants. The closure explicitly accounts for the influence of the Damkohler number 
and, of course, includes the mixing solution in the limit of zero Damkohler number. Several 
closures for the pressure-scalar gradients correlations are considered and the predicted results 
are compared with available experimental data in nonreacting and reacting turbulent shear 
flows. 


II Theoretical Background 

With the convention that the angle brackets, ( ) represent the ensemble mean value of a 
transport variable and the prime denotes its fluctuations from the mean, the nondimension- 
alized averaged equations for incompressible, isothermal turbulent reacting flows are: 


dxj 


= 0 , 


(i) 


a( «.) a{ u,){ Uj ) _ wq l M i a>.) 

dt dxj dxj (p) dxi Redijdxj' 

i,j = 1,2,3 (2) 
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( 3 ) 


d(Y a ) d^juj) dW® 1 0*(r o ) . , 

dt dxj dxj ScRedxjdxj ' ah 

a = A,B. 

Here u,-, P, p, Y a , Re, and Sc denote the i-th component of the velocity vector, the pressure, 
fluid density, mass fraction of species a, the Reynolds number, and the Schmidt number, 
respectively. { u > a ) represents the rate of chemical reaction = (<^b))'- 


<«.) = -Da{{Y A )(Y B ) + («) 

where Da is the Damkohler number. The algebraic formulation involves a two-equation 
scheme in which the Reynolds stresses and the scalar fluxes are expressed by nonlinear 
functions of the mean gradients and the time scales of the flow. The mechanical time scale 
is determined by the solution of transport equations for the kinetic energy of the turbulence 
(Jfc) = (ui u i )/ 2 an d for the turbulent dissipation: (e) = For boundary-free shear 

flow, these equations are: 9,27 


m , gwK) 

dt dxj 


■w, + V)^. 


+ 


1 d 2 {&) 

Re dijdxj 


_ I/MM) 

' • >’ dxj Re''dx j dx j 1, 


( 5 ) 


dt dxj 


_d_ 

dxj 


. , . 2 dp' du'- 

^ + lle^dxkdx,) 


+ 


1 a 2 (e) 


Re dxjdxj 

Cti {ky i j 'd Xj Gi2 (k) ’ 


( 6 ) 


with C t , = 1.4 and C t3 = 1.9. Treatment of the scalar variable requires the solution of 
additional transport equations for the reactants’ covariance (Y^Yp), and dissipations {tap) = 
’ScRe^Mo L ^') ■ F° r f° rrner we have: 
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d{Y'Yj) dKYfiWi) aWjYjYS) I 1 PPM) 

dt dxj dij ScRe dijdxj 

,,y^{Y.) _ 2 aYjdYj 

' J 13 dxj ScRe ' dij dxj 

+R*Z> + (<W- (?) 

By neglecting the third-order mass fraction correlations, the chemical source terms in the 
expanded form read (no summation on greek indexes) (u> a Yp) + (upY^) = —Da[[{Y^Y^) + 
{YpYa))(Yb) + + OtP'bJJO'j*) ]• Full resolution of the nonlinear interactions in the 

chemical scalar fields requires significant computational effort in practical applications. The 
neglect of the higher-order scalar fluctuations is justified by earlier results 28 but cannot be 
recommended for general applications. In such applications, the single-point probability den- 
sity function (pdf) or the joint pdf of the scalar variable provides the required information. 29 
The inclusion of the pdf is not attempted here. 

By using the first order term in the two-scale direct-interaction approximation, Yoshizawa 26 
develops a generic model for the scalar dissipation. An equivalent expression is obtained from 
Yoshizawa’s results (by making use of the dissipation time scales ratio r a = e(Y^)/(k)e a ). 
The equivalent form of this equation including the effects of chemical reactions is of the 
form: 30 


■WX) 


m) 

dxj 


d{ e ap) d{e a p)(uj) _ d . / . 1 d^[tgp) 

dt dxj dxj ' i 0,13 ScRe dxjdxj 



In the limiting case of pure mixing, the magnitudes of the constants are: C Vl = 1.7, = 1.4, 

C V3 = 2.0 and C Vi = 0.9. Similar values are suggested in Ref. 31 In this equation, the chemical 
source term source is of the form: 


S* = -Da [((£.) + {e„«))(Ji) + «erf + (e ol ,))(n)], a#/). (9) 
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Each of the reactants’ variances obey a similar transport equation with the chemical source: 


S aa = -2 Da ((e a )(Y x ) + (e aA ))(y a )) , (10) 

where A = B if a = A and vice versa. To complete the closure formulation, all the third- 
order transport terms are described by the gradient diffusion hypothesis. Denoting by 0 any 
of the second-order scalars, we have: 


K6) = -c.$ 


( 11 ) 


where C s is taken to be equal to 0.22 for all correlation type of quantities, whereas for the 
turbulent dissipations, C s = 0.18. Also, the molecular transport terms are neglected under 
the assumption of high Reynolds- Peclet numbers flow. 


Ill Explicit Algebraic Models 

An improved explicit algebraic Reynolds stress model (ARSM) has been derived by Taulbee 2 
from the modeled transport equation for the Reynolds stress. This model is based on the 
general linear pressure-strain closure given by Launder et al. 27 The improvement is due to an 
extended range of validity; the model is valid in both small and large mean strain fields and 
time scales of the turbulence. In two-dimensional flows, Pope 1 was first to give an explicit 
solution for the standard ARSM. The analogous improved ARSM as developed by Taulbee 2 
reads: 


a = —2 C M r + bigrcr 2 (^S — S^) + b 2 gT(SSl — flS)j . (12) 

Here, a = [a,j] = [(u^u'j) / (k) — 26ij/3] designates the anisotropic stress tensor, S = [5,j] = 
[d(ui)fdxj + d(uj)/dxi] and fl = = [d(ui) /dxj — d(uj) /dxi\ denote the mean flow strain 

rate tensor and mean flow rotation rate tensor, respectively. 6^ = [<$,^] = 1 for i = j = 1,2 
and 0 otherwise, r = (k)/{e) is the time scale of the turbulence and a = (SijSji) 1 ^ 2 is a 
strain rate tensor invariant. The parameters C M and g are given by: 
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9 = [Ct + - 2 + (2 - C tl )P/(e))-\ (13) 


r = 4^/15 

M 1 - Kbigrcr) 2 - 2(b 2 gruj) 2 ’ 

where P = —(k)aijSji is the production of turbulent kinetic energy, u> = (fl.jflj,) 1 / 2 is 
a rotation rate tensor invariant, and C\ t b\ and b 2 are constants from the pressure-strain 
correlation model (£>i = 0.086, b 2 = 0.377). Lasher & Taulbee (private communication) 
reffitted the expression for C\ proposed in Ref. 8 to six return-to-isotropy experiments and 
four homogeneous shear flow experiments obtaining for the linear model: 

Ci = 1.0 + F/18.0 exp(— 3.1/^/flej)(29.0 log(1.0 - 17.75(7/ + 7.1577/)) 


where F = 1 + 27777/8 + 977/4 is a parameter involving the second invariant 77 = 

and third invariant 77 = of the Reynolds stress anisotropy tensor, and Re\ = 

4(A:) 2 /(9(e)i/), the Reynolds number of the turbulence. 

A similar line of reasoning is made to obtain an algebraic closure for the velocity-scalar inter- 
actions. The transport equations governing these correlations are transformed into algebraic 
expressions by making two primary assumptions: (1) Existence of a “near-asymptotic” state, 
and (2) A negligible difference in the transport terms. The starting equations for the con- 
vective scalar fluxes are described by: 


MO . = gjKgTO + VVDMiy) 

dt dxj dxj 

+ V)( P 'S-{ { ^W + ^ K) W 

-Da((uX){yp) + {u'X){Y„) + JJJ)) 

_L \JL ( Y '?< + 

Re dxj \ a dxj ScdxjJ ScRe \ dxj dxj / 


+ 


(14) 


In the rhs of this equation, the following terms are identified: turbulent transport, pressure- 
scalar gradient correlation, production by the mean velocity and the mean scalar gradients, 
chemical reaction effects, molecular transport (assumed negligible at high Peclet numbers) 
and viscous dissipation. Based on the Poisson equation satisfied by the pressure fluctuations 
one can arguably split the pressure-scalar gradient correlation into two parts corresponding 
to so called rapid and slow terms. 8 The rapid term represents an inner product between the 
velocity gradient tensor and a third-order tensor, the last one subject to symmetry, continuity 
and normalization constraints. As suggested by Lumley, 8 since the slow pressure-scalar 
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gradient term and the viscous dissipation term are functions only of turbulent quantities, they 
can be incorporated into a single closure. The ensemble of the entire pressure-gradient term 
and viscous dissipation term enjoys a general relation encompassing some of the formulations 
proposed in precedent works. This is written consequently: 


*. = J_ 

” <P) 


JY'\ l+Sc/dvtNZy Ci g (e) 


dxi / ScRe \dxj dxj 


2 (k) 




+ 


d(ui) - fy.,. d(uj) . , t . d(xij) // y/v 

jYa) + + C3 ~d^ a ' AukYa) 

+Ct ^f teteW + a MK» + 


(15) 


The coefficients can take the following values: C la = 6.4, C\ = 0.5, c, =0, i = 2, 5; 18 
Ci a = 18.0(1.0 -(- 130.0/5cile) o - 25 (1.0 4- 12.5 /Re° A8 )~ 208 , c, = 0., i = 1,5, where Re t = 
4(A:) 2 /«c) I /); 32 or Cj = 4/5, c 2 = -1/5, c 3 = 1/10, c 4 = -3/10, c 5 = 1/5. 33 


To proceed, let us denote the 
flux) by: 


mechanical-chemical correlation coefficient (normalized scalar 



(16) 


where (\P) = (( Y£ ) + (Yg ))/ 2 is the turbulent scalar energy. The results of direct numerical 
simulations of homogeneous shear flow 32 suggest the existence of an asymptotic state for 
the normalized correlation coefficient ipi a , but not for the scalar flux itself. This observation 
justifies our first assumption. The second approximation is yet to be substantiated. The 
scalar flux equation can now be represented in terms of the correlation coefficient <fi a ( a ^ P), 
written symbolically: 


V^»’a 

IT 


D<p ia _ 
Dt 

(¥) dl V 

^ ( k ) dxj 


1 


Via 

«<*> 

ydxj 

2 ^ 


{, k ) dTf 

N (®) dx > 


Wafa) 

2* 


( — — 14- 

\fa) fa)) 


Via 
2 T 



4 P ia + + Sic, 


(17) 
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where the shorthand notations D/Dt, T?, Tf and Tj stand for convective transport, tur- 
bulent transport of scalar flux, scalar energy and kinetic energy, respectively. Moreover 
Po = -VWlMft)/*'/ + VjBd(Y B )/dx j) is the production of scalar energy, (c#) = 
(za) + ( c s) denotes the dissipation of scalar energy, S<s = {Cj^Ya) + {ubYq) is tbe chemical 
source term in the 'f’ equation, and the remaining quantities are the normalized production, 
pressure-gradient and the chemical source term: 


Pic = - 


(k) 2 d(Y a ) 

,j + 3^ dxj + fi u) 


(18) 


c 

^ ta — ^ Pta [(^1 “1" ^2)SijP jet + (Cx 

Zt 

+ (^3 + C\)d{jS jkpka ^S^jk^ijPka “ I " (^3 ^4)^ij^jkpka 

+ Csajk$lij<Pka + C^djkSjkPia] 


(19) 


$ ia — -Da{Vi Q (Y 0 ) + ipip(Ya) + liap\j (^))- 


( 20 ) 


Here 7 ia p = {u',Y^)/ yj(k)(<i>) 2 is the normalized covariance flux vector. Under the stated 
assumptions, the terms representing the convective transport and the turbulence diffusion 
difference are set to zero. This procedure leads to an algebraic system of equations for the 
two unknown vectors <p, a and ipip: 


<p a + D a A.(p a + B a <pp + C 0 = 0 
ifp + DpAipp + Bpip a + C p = 0 


where the coefficients are 


D a = 


2rft a 

TTw^hJXp)' 


Dp = 


2rhp 

1 + 2Dahp(Y a ) 


B a = Da(Y a )D a Bp = Da(Yp)Dp , 


with 



£lcr — 1 + (1 + 2c4)-^-y 


^(e») / P* . S<t \1 1 

<*) \{e*) i + (e*))\ ’ 


( 21 ) 


( 22 ) 

(23) 


( 24 ) 
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and the free vector terms read: 


C ia = 


2 rh a (k),_ ,2 , ,d(YJ 

i + 2Darh a (Y fi )\ (9) {a * + 3 ki) dx k 


2Th a DayJlj$) 

1 + 2DaTh a (Yp) 


IfiotP > 


( 25 ) 


„ 2r h, <*> . 2 C ^(lj) 

C '° - 1 + 2DaThjj(Y a ) \ <») ( “ + ) “ 3x t 


2rhpDayJ{^) 

1 + 2DaThp(Y Q ) 


7«dr0* 


( 26 ) 


Finally, the anisotropy of the turbulent diffusivity is ensured by the properties of the second 
order tensor A: 


= [(1 — ci — c 2 )S ik + (1 — cj + c 2 )ft,-jfc — (c 3 + c 4 )aijSjk 

— c$o.kjSji (c 3 c 4 ')cLijH,jk -t - c^cikjH jj] • (27) 

This tensor turns out to be traceless (An = 0) as a consequence of incompressibility and 
of the particular values taken by the constants c,-’s. Now, the solution of the system is 
conveniently represented in the matrix form: 

I V. = -M- 1 [(« + D 0 A)C. - B„ C„] 

( (pp = — M -1 [(6 + D a A)Cp — BpC a ] 

where M denotes the matrix [(1 - B a Bp)S + (D a + Dp) A + D a DpA 2 ]. The new expressions 
for the turbulent fluxes of reacting scalars exhibits two novel features. First, there is a direct 
influence of the Damkohler number Da, hence one can expect a higher flexibility of the 
model with respect to chemistry. Second, the strong coupling existing between the evolution 
of the reactants is reflected by the nonlinear dependence between the mean scalars and the 
covariance flux in the closure. 

To provide a computationally efficient algorithm, the matrix M is inverted analytically. This 
is achieved via the use of the Cayley-Hamilton theorem and yields a vectorial expansion 
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defining a natural basis for this problem: 


2 2 

Va = ^2 a nA n C a + ^ a( l A n C , / 5 . (29) 

n=0 n=0 V 

In the Appendix the inversion procedure via the Cayley-Hamilton theorem is outlined and 
the coefficients a n and a' n are listed. The final results provide an explicit solution for the 
scalar fluxes. In the limit Da -► oo, the coefficients a n become singular. In this limit, it 

is recommend to use the mixing solution Da = 0 for the transport of a Shvab-Zel’dovich 
variable. 34 


IV Illustrative Examples 

In this section, some sample results are presented of numerical simulations based on the mod- 
els presented above. There are two primary reasons for conducting these simulations: (1) 
Model assessments via comparisons with laboratory data, (2) Demonstration of the model 
capabilities in comparison to traditional closures based on the linear gradient diffusion ap- 
proximation. The flow configurations considered are those for which abundant laboratory 
data are available: two-dimensional planar mixing layers under both nonreacting and non- 
premixed reacting conditions, planar jets and axisymmetric jets. The mean flow motion in 
all of these shear flows is assumed two-dimensional. The space coordinates are identified 
by x = [x,y], x is the streamwise coordinate denoting the direction of principal evolution 
of the flow, and y represents the cross-stream direction. The velocity field is identified by 
U = [«,»]. In nonreacting flow simulations the mass fraction of one conserved species, V A 
is considered. In the mixing layer configuration, V A = 1,0 at the low-speed ((u) = u L ), 
and at the high-speed ((u) = u H ) streams, respectively. The magnitudes of r = u H /u L 
and/or Au = u H - u L are set according to each of the experiments considered. In the jet 
configurations, Y A = 1 is issued from the jet into a surrounding of V A = 0. In the reactive 
mixing layer, the two species A and B are introduced into the low- and high-sped sides. 
These species are assumed thermodynamically identical and there is no trace of one of these 
species at the feed of the other one, i.e. complete initial segregation. The mass fractions of 
the two species at their respective feeds are set equal to one. In accord with the reacting 
flow experiments, the heat generated by the reaction is assumed negligible. 

The transport equations governing the velocity and the scalar fields are of parabolic type 
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with the boundary-layer approximation. The numerical algorithm is based on a first-order 
upwind differencing for the convection terms and a second-order central differencing scheme 
for all the other terms. Due to the anisotropic character of the algebraic closures, it is possible 
to evaluate all the components of the Reynolds stress tensor and the scalar flux vectors. The 
Reynolds stress tensor is predicted by the explicit ARSM solution as developed by Taulbee. 2 
The scalar flux vectors are computed with the solution given by Eq. (28). In this equation, 
several different closures as given in Refs. 18,32,33 are used for the pressure-scalar gradient 
correlation. 

In the evaluation of the Reynolds stress tensor and the scalar flux vectors, the terms appear- 
ing as model coefficients ( e.g . P/(e) in Eq. (24)) are treated as known quantities. Therefore 
an iterative numerical solution procedure is required. To insure a faster convergence and to 
avoid numerical instabilities, the initial profiles are set close to those obtained by the simi- 
larity solution. The implementation of the boundary conditions is similar to that in many 
previous simulations of parabolic shear flows (e.g. Ref. 28 ). In the results presented below, 
the spatial coordinates are presented by 77 for hydrodynamic and 77* for the scalar variables. 
In the mixing layer, 77 = y-y ff g 0 °' 5 ^ , 77* = y ~ y ^^ = °' 5 ^ where u = and x 0 is the virtual 

origin. In the jet configurations 77 = and 77* = 5 j . In these jets, the subscript 

CL denotes values at the centerline (i.e. y = 0). In the reacting mixing layers, the corre- 
sponding profiles of Y a under no chemistry is used in the normalization of the cross-stream 
coordinate. 

In Figs. 1-4, results are presented of nonreacting mixing layer simulations and are com- 
pared with data obtained in several laboratory experiments. 35 ^ 11 In these and all the figures 
presented below, the transverse variations of statistical variables are presented. In this par- 
ticular flow, the model with the coefficients of Ref. 32 yields a constant C\ a which is too 
small for the numerical stability; the closure of Ref. 33 provides the best overall agreement. 
Figures 1-4 indicate that the agreement between the model predictions and laboratory data 
is generally good for the mean values of the streamwise velocity and the mass fraction of 
a conserved scalar. The exception is the behavior near the high-speed stream of Masutani 
and Bowman’s 35 layer where the experimental results are underpredicted by the model. A 
similar trend has been observed in previous simulations based on a linear gradient diffusion 
model. 28 A larger skirt of the mean scalar profile is believed to be due dominant influence of 
large scale structures in this particular experiment. 42,43 The model yields a better agreement 
with the results of Batt’s 36 experiment, in which the flow is less dominated by large scale 
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structures. For the second order moments, the scatter in reported data does not allow a 
precise comparison. In general, the trends are in accord with data and the agreement can 
be potentially improved by modifying the magnitudes of turbulent scales at the boundaries. 
Associated with the plateau in the measured mean scalar values of Masutani and Bowman, 35 
there is a double hump in the experimental scalar variance (Fig. 4(a)). The locations of these 
humps coincide with regions of localized large mean scalar gradients. This behavior is cor- 
roborated by direct simulations 42,43 of two-dimensional, weak turbulent flows but cannot be 
predicted by the model and is not observed in large Reynolds number experiments. 36 

The anisotropy of the Reynolds stress tensor and the scalar flux vector as predicted by 
our model, allow a direct comparison of the predicted fluxes with data. Again, the general 
agreement with laboratory data as witnessed in Figs. 5-7 is satisfactory. Moreover, with these 
results it is possible to perform posteriori appraisal of the closures based on conventional 
linear gradient diffusion hypothesis. For example, the parameters and Sc t as given by: 

/ / a d(u) „ (k 2 ) , x 

(uv) = (30) 

v'Y') = — QQliI 

a) Sc t dy ’ 

can be directly evaluated. The explicit algebraic relation for C M is given by Eq. (13); the 
relation for the turbulent Schmidt number obtained from the model with coefficients of Ref. 33 
is: 



2g 1 - [ a ?2 + 3 (fl22 + |) (2 — an — 2 a 22 )] 

15/,. [j + _ §)( S r^)’] (l - (<■» + 1) ' 


Figures 8-9 show the streamwise variations of these parameters. A somewhat similar qualita- 
tive behavior is observed in the jet flows. These results can be compared with = 0.09 and 
Set = 0.7 typically employed in the linear gradient diffusion approximations. 44,6 Also, the 
ratio of the scalar to velocity time scales (r 0 ) as shown in Fig. 10 indicates that an approx- 
imate constant value can be used for the central region of the layer. This is in accord with 
the results in Refs. 45,46 However, there are large peaks near the free-streams. The amplitude 
of the peaks can be somewhat controlled by modifications of the boundary conditions. An 
exact specification of these conditions requires inputs from the laboratory measurements. 
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The results for the reacting mixing layer are presented in Figs. 11-12, where comparisons are 
made with experimental results. 35 Figure 11 indicates that the predicted mean mass fraction 
of reactant A is in accord with data. Because of the consumption of the reactant by the 
chemical reaction, the discrepancy near the free-stream as noted in Fig. 3 is not observed here. 
Naturally, the second hump in the experimental profiles of the reactant’s variance vanishes 
as the species is consumed. The variance amplitude is, nevertheless, still underpredicted by 
the model as shown in Fig. 12. The cross-stream variation of the reactants’ covariance (the 
unmixedness) as predicted by the model is shown in Fig. 13. In accord with the physics 
of turbulent flows with segregated reactants, the unmixedness is negative throughout the 
layer. The same is true in the limit of zero chemistry (mixing only) as shown in this figure. 
However, in this case the amplitude of the unmixedness is slightly larger. This is not in 
agreement with Toor’s 47,48 hypothesis which suggests an independency of unmixedness to 
the Damkohler number, but is consistent with the more rigorous theories of nonpremixed 
turbulent reacting systems 49-52 which indicate dependence on the Damkohler number. The 
results in Figs. 8-10 also suggest that chemistry induces a profound influence on the amplitude 
of the model constants. 

The measurement results pertaining to planar jet hydrodynamics as reported in Refs. 53-55 are 
compared with the model predictions in Figs. 14 through 20. The agreement is reasonable 
for the mean streamwise velocity and also the components of the Reynolds stress tensor. The 
same level of agreement is witnessed in the profiles of the conserved scalar mean and variance 
as compared with the experimental results of Refs. 56-60 In this case, it seems that the model 
with the coefficients of Ref. 33 provides the closest prediction of the scalar mean and variance 
values. However, the predicted values of the amplitude of scalar-velocity correlation is larger 
than some of the experimentally reported data. 

The comparative assessment of the models for the prediction of axisymmetric jet flows is 
provided in Figs. 21-27 where the experimental data in Refs. 61-63 are used for hydrodynamics 
variables and those in Refs. 64-66 for scalar variables. Again, all the mean values are reasonably 
well-predicted. The same is true for the Reynolds stresses and the scalar variance as the 
model results are within the scatter of the data. In this configuration also the experimental 
results for the scalar fluxes are overpredicted. However, the experiments of Ref. 64 are not 
conducted in the self-preserving regions of the jet. Therefore a definite verdict cannot be 
reached without comparisons with further data. 

From the preceding comparisons, it can be concluded that the algebraic model developed 
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here provides an effective means of predicting the second-order moments in reacting turbulent 
flows. Because of their anisotropic feature, these algebraic schemes are more flexible than 
the conventional linear gradient diffusion schemes. The explicit nature of the relations, 
facilitated by the Cayley-Hamilton theorem, is particularly convenient for simulating complex 
flows. Additional improvement of these models is possible by fine-tuning of the pressure- 
gradient correlation models, and also the higher order moments of scalar-scalar fluctuations 
in reacting flows. Further extension of the model is recommended to account for the effects 
of exothermicity in nonequilibrium chemically reacting systems. 


Appendix 

The procedure leading to explicit solutions for the scalar-flux vector, as governed by Eq. (21) 
is described here. 

Consider an arbitrary three-dimensional second-order tensor Q = [Q.j] and the correspond- 
ing Kronecker tensor 6 = [<5 tJ ]. According to the Cayley-Hamilton theorem, this matrix 
satisfies its own characteristic polynomial: 

Q 3 - IqQ 2 + IIqQ - IIIqS = 0 (33) 

where I Q = {Q} = <?,,, IIq = |[{Q} 2 - {Q 2 }] = \[Q«Q» ~ QijQji], IIJq = |[{Q} 3 - 
3{QHQ 2 } + 2{Q 3 }] = liQaQjjQkk - 3 QaQjkQkj + 2 QijQjkQki] are the three tensorial 
invariants. Multiplying the characteristic polynomial with Q -1 and solving for the inverse 
we obtain: 

Q-' = ~(Q , ~ hQ + HqS). (34) 

This relation can be used now to find explicit solutions to the problem considered here. For 
example, for the case with Da = 0 (pur mixing) we can write: 


<p a — —(6 + G) 2 C, 


(35) 


where G = D a A. Hence: 


(6 + G) -1 = 


1 

TTi 5+G 


[G 2 + (2 — I s+ g)G + (1 — Is+g + Ils+an 


(36) 


14 


It is easy to show that I s+ g = Ig + {<$}, Ih+G = 2/g + Ug + {£}> Uh+G = Ig + Ug + 
IIIg + 3 • Therefore the normalized scalar flux vector takes the form: 

<f a = oqC + aiGC + a^G^C (37) 


with the coefficients: 


1 + Ig + 7 /q 
1 + Ig + IIg 4 * IIIg 
1 + Ig 

1 + Ig + Ug + IIIg 
1 

1+Ig + IIg + IIIg' 


(38) 

(39) 

(40) 


In a homogeneous solenoidal field the pressure-gradient correlations are described by traceless 
rapid parts. This further translates into {A} = 0, and thus: 


1-I{G’} + I{G3} 


(41) 




a 2 


1 

i-H G2 } + H G3 } 


i 

l-i{G 2 } + H G3 }' 


(42) 

(43) 


The reacting case is somewhat more complex. Nevertheless, by following the same procedure 
explicit solutions are obtained: 


<, p a — uoG q + a' 0 Cp + ai AC 0 + a\ACp + a 2 A?C a + a' 2 A 2 Cp (44) 

V>fi = b 0 C a + b' 0 C p + 6,AC q + b[ACp + b 2 A 2 C 0 + b 2 A 2 Cp, (45) 

with the coefficients: 

Fp(F a + Dl *£1) + B a Bp [Up-(D p F a - ED a )Dp - E(E + ^p-D a Dp)} 

°° (1 -B a B 0 ){F a F 0 -E 2 B a B 0 ) ’ (46) 


, F Q (Fp + D}&1) + D a ^(D a F 0 - EDp ) - EB a Bp(E + ^D a D 0 ) 

(1 -B a Bp)(F a F p -E 2 B a Bp) 


(47) 


BgBp^ [E{Dp — Do) + FgDp\ — DgFp 
D a (F a Fp - E 2 B a Bp) 


(48) 
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(49) 


a\ - 


D a 


a 2 = 


CZo — 


with the shorthand notations: 


F a = (1 — B a Bp)(D, 


Fp = (1 — B a Bp)(Dp 


— DpF a 

— E(Dp — D 0 

F a F p 

- E 2 B a Bp 

D a Fp — 

EDpB a Bp 

F a Fp- 

E 2 B a Bp 

n D a Fp — EDp 


- E 2 B a Bp ’ 

. {A 2 } 

| 

£ 

0 2 

D a Dp 

, (A’> 

1 B a Bp 

0 2 

Dp D a 


{A’} 


) 2 

“ 3 


)-■»; 


, {A 3 } 


E = (5- + i)(l - B„B e ) - DvDp^p-. 

The coefficients fcj are obtained from the a,’s through the permutations a 
ao — ► 6q, aj, — »■ 60 , ai — » 6^, a'j — > 6 j, 02 — * ^2 a 2 b 2 . 


(50) 

(51) 


(52) 


(53) 


(54) 
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Figure Captions 

Figure 1. Cross-stream variation of u for the mixing layer. 

Figure 2. Cross-stream variation of (u ,2 ) 1/2 for the mixing layer. 

Figure 3. Cross-stream variation of (1*) for the nonreacting mixing layer. 
Figure 4. Cross-stream variation of {Y£) l/2 for the nonreacting mixing layer. 
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Figure 5. Cross-stream variation of (u'v')/Au x 100 for the mixing layer. 

Figure 6. Cross-stream variation of {v' 2 ) 1 ! 2 j Au for the mixing layer. 

Figure 7. Cross-stream variation of (v'Y A ) / Au for the nonreacting mixing layer. 

Figure 8. Cross-stream variation of C ^ for the mixing layer. 

Figure 9. Cross-stream variation of Sc t for the nonreacting and reacting mixing layers. 
Figure 10. Cross-stream variation of ta for the nonreacting and reacting mixing layers. 
Figure 11. Cross-stream variation of (Ya) for the reacting mixing layer. 

Figure 12. Cross-stream variation of { Y 2)'^ for the reacting mixing layer. 

Figure 13. Cross-stream variation of (Y A Yg) for the nonreacting and reacting mixing layer. 
Figure 14. Cross-stream variation of ( u)/(u)cl for the planar jet. 

Figure 15. Cross-stream variation of («’ 2 )/(u) 2 cl for the planar jet. 

Figure 16. Cross-stream variation of (uV)/( u )cl for the planar jet. 

Figure 17. Cross-stream variation of (t/ 2 ) /(u)q L for the planar jet. 

Figure 18. Cross-stream variation of (Ya)/(Ya)cl for the planar jet. 

Figure 19. Cross-stream variation of {Y A ) 1 ^ 2 /{Ya)cl for the planar jet. 

Figure 20. Cross-stream variation of {v'Y a ) / (u)cl(Ya)cl for the planar jet. 

Figure 21. Cross-stream variation of (u)/(u)cl for the round jet. 

Figure 22. Cross-stream variation of (u 13 ) I (v)c L for the round jet. 

Figure 23. Cross-stream variation of (uV)/( u )cl for the round jet. 

Figure 24. Cross-stream variation of (v 12 )/ (u)q L for the the round jet. 

Figure 25. Cross-stream variation of (Y a )/(Ya)cl for the round jet. 

Figure 26. Cross-stream variation of (Y^) 1 ^ 2 /(Ya)cl for the round jet. 

Figure 27. Cross-stream variation of (v'Y a )/ (u)cl(Ya)cl for the round jet. 
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Abstract 

Experimental data of Bilger et at. (1991) pertaining to the compositional struc- 
ture of a turbulent reactive scalar mixing layer are reproduced by a mathematical- 
computational procedure utilizing the Pearson family (PF) of univariate and mul- 
tivariate scalar probability density functions (pdfs). Aided by comparisons against 
these data and some additional data generated here by direct numerical simulation 
(DNS) of a spatially developing reacting mixing layer, an appraisal is made of the 
applicability and the extent of validity of the PF for the statistical description of the 
reactant fluctuations. In accord with the experiment, a chemical reaction of the type 
A + B — ► Products is simulated in isothermal, incompressible flows. A wide range of 
the Damkohler number is considered including both frozen and equilibrium chemistry 
flows. The comparison of the predicted results with laboratory data indicates that the 
PF generated pdfs are very convenient, in the absence of better alternatives, for mod- 
eling the influence of turbulence on the mean reactant conversion rate. In particular, 
the Dirichlet frequency parameterized with the “scalar-energy” provides the most rea- 
sonable means of portraying the multivariate scalar pdf. A more detailed comparative 
assessment of the model predictions against the DNS data confirms the relative supe- 
riority of the Dirichlet pdf even though it is mathematically shown that this frequency 
cannot be justified for modeling of equilibrium flows. 

With the use of the PF generated pdf, the autospectral density function and the 
cross-spectral density function of the reacting scalars under equilibrium chemistry are 
related to the frequency spectrum of the mixture fraction. This relation is very conve- 
nient and is favored over previous models for predicting the spectral characteristics of 
reacting scalars in the central region of the laboratory mixing layer and in any other 
homogeneous flow configuration which yields an asymptotic gaussian pdf for the mix- 
ture fraction. The correction to these spectral relations for asymptotic exponential 
pdfs is provided, and the influence of the Damkohler number on both the auto- and 
the cross-spectral density functions is assessed by the analysis of the DNS data. 
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1 Introduction 


In a recent article, Bilger et al. (1991) (hereinafter referred to as BSK) report the results 
of extensive laboratory measurements on the compositional structure of a turbulent reactive 
scalar mixing layer. The objective of these measurements is to portray the behavior of the 
reacting scalar field in both “single-point spatial” and “two-time” statistical levels. In the 
format presented, the data are very useful for assessing the influence of scalar fluctuations 
on the rate of mean reactant conversion and also for examining the spectral characteristics 
of reactive scalars in a chemical reaction of the type A + B -* Products. A wide range of 
chemistry parameters is considered including both frozen and equilibrium flows. The extent 
of the measured data together with the relative simplicity of the flow configuration provide 
an excellent setting for appraising the performance of turbulence closures in accounting 
for the role of scalar fluctuations (Komori et al., 1991), and for predicting the spectra of 
reacting scalars. The conclusion of this work advocates the need for mathematical models for 
depicting the stochastic and the spectral character of reacting scalars for general applications. 

Mathematical modeling of scalar fluctuations in stochastic treatment of turbulent react- 
ing flows has been the subject of broad investigations since the pioneering work of Toor 
(1962). One approach which has proven particularly useful is based on the probability den- 
sity function (pdf) or joint pdf of scalar quantities (Libby and Williams, 1980; O’Brien, 1980; 
Bilger, 1980; Pope, 1985; Givi, 1989; Kollmann, 1990; Pope, 1990). This approach offers the 
advantage that all the statistical information pertaining to the scalar field is embedded within 
the pdf. Therefore, once the pdf is known the effects of scalar fluctuations are easily deter- 
mined. Because of their capabilities, pdf methods have been widely utilized for a variety of 
reacting systems (see Libby and Williams (1994); Pope (1994) for recent reviews). A system- 
atic approach for determining the pdf is by means of solving the transport equation governing 
its evolution (Lundgren, 1967; Lundgren, 1969; Pope, 1985). In this equation the effects of 
chemical reactions appear in a closed form, but modeling is needed to account for the trans- 
port of the pdf in the composition domain of the random variables. This transport describes 
the role of the molecular action. In addition, there are extra dimensions associated with the 
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composition domain which must be treated. These problems have constituted a stumbling 
block in utilizing pdf methods in practical applications; developments of turbulence clo- 
sures and numerical schemes which can effectively deal with these difficulties have been the 
subject of numerous investigations within the past two decades (Libby and Williams, 1980; 
Libby and Williams, 1994). 

An alternative approach in pdf modeling is based on “field-parameterization” methods. In 
these methods the pdf is not determined by solving a transport equation. Rather, its shape 
is “ assumed ” in terms of the low order moments of the random variable(s). Obviously, this 
method is ad hoc but it offers more flexibility than the first approach. Presently, the use 
of parameterized methods is justified in cases where there is strong evidence that the pdf 
adopts a particular distribution (Bilger, 1980; Pope, 1994). 

Between these two approaches, obviously the first is preferable if an appropriate closure 
is available to account for the effects of molecular action. In its application in turbulent 
combustion, traditionally the family of models based on the coalescence/dispersion (C/D) 
closures (Curl, 1963; Janicka et al , 1979; Pope, 1982; Kosaly and Givi, 1987; Norris and 
Pope, 1991), or linear mean square estimation methods (Dopazo and O’Brien, 1976; O’Brien, 
1980) have been employed. These closures are plausible from a computational standpoint 
and can be effectively simulated via Monte Carlo numerical methods (Pope, 1981). However, 
there are several drawbacks associated with these closures which restrict their use for reliable 
predictions (Pope, 1982; Kosaly and Givi, 1987). Some of these drawbacks are overcome in 
the newly developed Amplitude Mapping Closure (AMC) (Kraichnan, 1989; Chen et a/., 
1989). This has been established in a number of recent validation assessments of the AMC 
by means of comparison of its predicted results with those of DNS (Pope, 1991; Madnia et 
al., 1992; Jiang et al., 1992; Frankel et al., 1993), and experimental (Frankel et al., 1992) 
data. 

Despite its demonstrated properties, there are some deficiencies associated with the AMC 
which require further investigation. These are discussed in detail by Miller et al. (1993); 
the most serious of these are: (1) the “single-point” nature of the closure, (2) the difficulties 
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associated with its numerical implementation especially in multivariate statistical analyses, 
and (3) its inability to account for the migration of scalar bounds as mixing proceeds. The 
first problem is shared with C/D models and indicates the deficiency of the approach in 
accounting for the variation of turbulent length and time scales. The other problems are 
exclusive to the AMC and can cause difficulties in its applications. 

Considering the current state of progress in pdf modeling, it can be cautiously argued that 
parameterized pdf methods are somewhat more “feasible” than the transport equation ap- 
proach for practical applications. This is not to suggest the superiority of assumed meth- 
ods. Rather, it is to encourage further research on the first scheme before it can be im- 
plemented routinely. In this regard, Miller et al. (1993) and Frankel et al. (1993) have 
conducted detailed investigations based on the two approaches. The general conclusion 
drawn from these studies is that in cases where the AMC has proven useful, other ap- 
proaches based on parameterized pdfs are equally effective. In the cases considered by 
Miller et al. (1993), it is shown that in simple flows where the AMC can be employed, 
the family of pdfs based on the Johnson Edgeworth Translation (JET) (Johnson, 1949b; 
Edgeworth, 1907) can also be used. In fact, for the simple problem of binary mixing in 
isotropic turbulence - a standard test case - the solution generated by the AMC (Pope, 1991; 
Gao, 1991) can be viewed as a member of the JET family. Furthermore, due to established 
similarities of the JET with simpler distributions belonging to the Pearson Family (PF) 
of pdfs (Pearson, 1895), it can be argued that the PF can also be considered as a viable 
alternative (Frankel et al., 1993). 

There is a long history of the use of the PF pdfs in turbulent combustion, t.g. Rhodes 
(1975); Jones and Priddin (1978); Lockwood and Moneib (1980); Peters (1984); Janicka 
and Peters (1982); for recent reviews see Williams (1985); Givi (1989); Priddin (1991). In 
most applications to date, this family has been used in the form of the Beta density of 
the first kind (Pearson Type I and Type II distributions). This is due to the flexibility 
of this density in portraying bimodal distributions. The properties of this density have 
been examined in detail by Madnia et al. (1992); Miller et al. (1993) and Frankel et al. 
(1993). According to these studies there are some similarities between the PF and the AMC, 
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as well as some distinct differences. As indicated before, both these methods are utilized 
in the context of a single-point closure. Therefore, in both cases the magnitudes of the 
moments up to second order must be provided externally. Also, both methods suffer from an 
inability to account for the shrinking bounds of the scalar domain as mixing proceeds. This is 
manifested by the failure of both closures in producing a correct evolution for the conditional 
statistics of the scalar; namely, the conditional expected dissipation and the conditional 
expected diffusion. This can be troublesome and may produce significant errors, especially 
in modeling of non-equilibrium flames. In modeling of equilibrium reacting homogeneous 
flows, both closures are satisfactory regardless of the equivalence ratio (Madnia et a/., 1992; 
Frankel et al . , 1993). However, the actual implementation of the AMC is very difficult, 
if not impossible, for applications in non-homogeneous flows regardless of the chemistry 
model. In such systems the mapping transfer function must be evaluated numerically, and 
in non-equilibrium flows the multivariate form of the pdf must be considered. These require 
serious investigations before they can be implemented routinely (Pope, 1991). In these cases 
the application of the PF is much more straightforward but obviously cannot be justified 
rigorously. The corresponding multivariate pdf is the Dirichlet frequency (Naxumi, 1923; 
Johnson, 1987; Johnson and Kotz, 1972; Wilks, 1962), and its mathematical implementation 
is straightforward. 

Description of the spectral characteristics of turbulent flows has also been a challenging issue 
since the early contribution of Taylor (1938). This is again due to the system non-linearities; 
and development of closures to account for the influence of convection in the spectral trans- 
port equations continues to be an active area of research (Stanisic, 1988; Lesieur, 1990; 
Orszag, 1977). Naturally, the degree of difficulty is escalated when modeling of the spectral 
density function of reactive scalars is attempted (Corrsin, 1961). In such modeling, in addi- 
tion to the influence of convective fluctuations the role of the chemical effects must also be 
considered. The challenge in developing reliable spectral closures for this purpose has also 
been long- recognized (O’Brien, 1960; Corrsin, 1981). Recent efforts are devoted on develop- 
ing closures by which the spectral transport equations are coupled with the single-point pdf 
evolution equation (O’Brien, 1985; Jiang, 1992; Frankel et al ., 1992). In these approaches, 
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the spectral closures include two-point statistical information, and the pdf provides the clo- 
sure to account for the influence of chemistry. 

The objective of this work is to develop mathematical turbulence closures for modeling of 
the fluctuations and the frequency spectra of reactive scalars in the context considered in 
the BSK experiments (and also in Saetran et al. (1989)). The goal is to analyze the influ- 
ence of the scalar fluctuations on the rate of mean reactant conversion and to determine the 
auto- and the cross-spectral density functions of the reactant fluctuations. In accord with 
the experiments, a homogeneous mixing layer is considered and the models are appraised by 
detailed comparisons against these laboratory data. Comparisons are also made with addi- 
tional data generated here by Direct Numerical Simulation (DNS) of a spatially developing 
planar mixing layer. The hope is to provide simple working relations for general applications. 


2 Flow Configurations 

A scalar mixing layer is formed when a traditional “grid-generated turbulent flow” is aug- 
mented with a scalar gradient (Libby, 1975; LaRue and Libby, 1981; LaRue et al ., 1981; 
Saetran et al., 1989). This can be facilitated by either a temperature gradient, i.e. by heat- 
ing half of the grids (Libby, 1975), or by introducing two different chemical species (BSK). 
The schematic of the homogeneous reacting scalar mixing layer is shown in Fig. 1 (a). This 
configuration consists of a unform mean streamwise velocity carrier inert gas contaminate 
by two reacting scalars, A and B. In the BSK experiments, doping gases are injected and 
are uniformly mixed in a carrier air stream. For the passive scalar experiments, nitric ox- 
ide (NO) is injected into one stream; in the reacting layer experiments, NO and O 3 are 
introduced in separate feeds into the air. The concentrations of the reactive species axe 
sufficiently small (of order of parts per million), allowing the assumption of constant density 
and isothermal flow. 

In the spatially developing mixing layer, in addition to scalar gradients there is also a cross- 
stream variation of the mean velocities. The schematic diagram of this configuration is 
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shown in Fig. 1(b). Two co-flowing streams traveling at different velocities are merged at 
the trailing edge of a partition plate. The reactants A and B are introduced into the high- 
and the low-speed streams, respectively. In this flow, in addition to small scale turbulent 
transport, the formation of large scale coherent vortical structures also plays a significant 
role on the rate of chemical reactant conversion. Therefore, the spatial inhomogeneity of the 
flow in the transverse direction makes the statistical analysis more challenging in comparison 
to that in homogeneous flows. 

There is a wealth of recent experimental data available pertaining to transport of reacting 
chemical species in spatially developing mixing layer, e.g. Koochesfahani and Dimotakis 
(1986); Masutani and Bowman (1986). However, the statistical results reported are not 
as detailed as those furnished by BSK. Therefore, alternatively, data axe produced here by 
means of direct numerical simulation. DNS has proven very useful in capturing many fea- 
tures of turbulent transport (Givi, 1994) and provides a powerful complement to laboratory 
measurements for the analyses of turbulent reacting flows. The configuration considered 
here is similar to that treated by Givi and Jou (1988); McMurtry and Givi (1992), but the 
statistical analyses are conducted in a format similar to those in the BSK experiment. 


3 Formulation, Modeling and the Computational Pro- 
cedure 

3.1 Modeling of the Reactant Conversion 

In both flows, all the chemical species are assumed thermodynamically identical and the 
fluid is assumed to be calorically perfect. The value of the molecular Schmidt number is 
set to unity and it is assumed that the two reactants are completely segregated at the 
inflow. With unity normalized mass fractions of the two species at their respective feeds, a 
mixture fraction denoted by / is defined (Williams, 1985; Bilger, 1980; Bilger, 1989) to have 
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normalized values of 1 and 0 in the streams containing A and B respectively. Under chemical 
equilibrium, the statistics of the reacting field are directly related to those of the mixture 
fraction. Therefore, univariate statistical analyses are sufficient. In non-equilibrium flows, 
the joint statistics of the reacting variables are required. The coordinates x and y denote, 
respectively, the streamwise and the cross-stream direction in both configurations (Fig. 1). 
Due to the single-point nature of the statistical treatment, all the moments up to second 
order serve as “inputs” for the pdf parameterization. In modeling of the homogeneous 
flow, the procedure developed by Libby (1975) is followed. Turbulent convective fluxes 
are modeled by the gradient diffusion closure with a constant turbulent viscosity v. The 
“self-similar” behavior of the layer is predicted within the domain of the transformed cross- 
stream direction, rj = YfyU*X , where Y = y/M, X = x/M ( M is the mesh spacing) and 
v * = vjUM (U is the streamwise velocity, and tilde denotes ensemble-averaging). With 
this framework, the modeled transport equations governing the ensemble mean value of the 
mixture fraction, /, and its variance, f' 2 (prime denotes fluctuations from the mean value) 
are given, in order, by: 


£l 

dri 2 



(i) 


<pr v ip 

1 dr} 2 2 dr) 


+ 2 



-C 2 f' 2 = 0, 


( 2 ) 


where Ci, C 2 are empirical constants. In the limit of infinitely fast chemistry, the pdf of the 
mixture fraction is sufficient to determine ail the statistical information pertaining to mass 
fractions of the species. The most convenient of the PF is the Pearson (1895) Type I in the 
form of the Beta density of the first kind: 


PfW = 


1 




( 3 ) 


Here, B(a,(3) is the Beta function: B(a,f3) = r(a)r(/?)/r(a + /3), T is the Gamma function, 
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and the parameters a,j3 are uniquely determined by / and f' 2 (Abramowitz and Stegun, 
1972). With this density, all the moments of the mass fractions of the reactants are expressed 
in a closed analytical form: 
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B(a,0) 

1 

WJ) 


1 1 

J (/• - /)”/“-'(! - df 

ft 

(1 £ (")(-/•)'/ - /)«-> df, 

ft 


YaYb = 0, (4) 

provided that n,m ^ 0. Here, YJ indicates the mass fraction of i-th species, and f, denotes 
the stoichiometric value of the mixture fraction. For unity normalized mass fractions at 
free streams, /, = 1/2. Analytic forms of these integrals are easily attainable and can be 
expressed in terms of the incomplete Beta density function (Madnia et al., 1992): 

T '- (a ' 0) =W)l° (5) 

In non-equilibrium flows, the transport equations for the first two moments (and/or cross 
moments) are coupled with the pdf. The transport equation for the ensemble-average of the 
mass fractions of the reactants, say species A is given by: 


<PYa r,jy A 

dr} 2 2 dr} 


Xu A =0 


( 6 ) 


where Cj a denotes the reaction rate of species A. The normalized form of this rate is given 
by Co a = DoYaYb ; Da is the Damkohler number (BSK). The closure for the evaluation of 
the mean reaction rate is provided by the pdf. Here, the Dirichlet frequency (Johnson, 1987) 
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is proposed for this purpose: 
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where V’i > 0, j = l,n; 51 V’j < 1; and £* > 0, fc = l,n + 1. The n+1 parameters (^’s) 

i=i 

can be evaluated if the values of an equal number of moments are supplied. Therefore the 
knowledge of at least one second order moment is necessary, a recognizable feature of single- 
point closures. It is straightforward to show that the product moments of order r = £ r j 

i - 1 

with respect to the origin are: 


/4,r, -r„ = J , V’n )dtMV>2 • • • dtp n 

Sn 
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where S n is the support of ■P(V’) in the composition space. With this, the ensemble mean 
value of the mass fraction of the z-th species is: 


Ji= (i 


71 + 1 1 

E (k 

k = 1 


( 9 ) 


and the r, + r_,th order correlation (product moment with respect to the mean values) between 
species rpi and xpj is: 


^TiTj 
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Therefore, the variance of species i and the covariance of species i and j are, in order: 


i>? = 


n+1 

fi(Efc-fc) 



n+1 n+1 

(EWU6 + 1) 

jt=i *=i 


; Wj = 




n+l n+1 

( E C*) 2 ( E 6 + l) 

k = 1 k=l 


(11) 


A particularly pleasing feature of the Dirichlet density is the condition J2 xj>j < 1. This 

j=i 

constraint is in accord with the law of conservation of species. However, the parameterization 
of the pdf is somewhat “overdetermined” in that it is not clear which combination of the 
moments is to be used to specify the parameters of the pdf (values of £,-). In the reaction 
scheme considered by BSK, n = 2; thus the parameterization requires the knowledge of 
three moments. The mean values of the two reactants, Ya,Yb are obvious choices; but the 
specification of the third moment is not clear. Here, we consider two options: (1) the scalar- 
energy Q = Y'a + Yg, and (2) the covariance (unmixedness) C = Y^Yg. These schemes are 
enacted with the relations: 


where: 


(i = Y a A; h = Y B A; 6 = (l-^-W; A = 6 + 6 + 6, 

x _ Ya(1-W + W-Yb) x 
Y? + Yg 


( 12 ) 


(13) 


with the first option, and 


A = 


YaY b 

Y'aY'b 


+ 1 ) 


(14) 


for the second one. This overdetermination is not encountered with the use of the multi- 
variate gaussian pdf, typically used in combustion modeling e.g. Bockhorn (1988a); Bock- 
horn (1988b). However this frequency has serious drawbacks in requiring an infinite support 
(Leemis, 1986) and not satisfying the law of conservation of species. The transport equations 
for the second order moments (and cross moments) are of a similar nature. For example, for 


II 



the unmixedness: 


<PYXYj y JYjXk I JMaWb 

1 dr) 2 2 drj dr) dr} 


- C 2 YXYi £ + X{Cj a YX + ubYX) = 0 . 


(15) 


The only additional modeling requirement is the specification of the empirical constants 
Ci and Ci- Here, the constants are determined in such a way as to yield the best agree- 
ment with the laboratory data for the first two moments, as will be discussed in Section 
4. The resulting system of nonlinear coupled differential equations is solved by means of 
the Newton- Raphson iterative scheme. The derivatives are approximated via a second order 
accurate finite-difference scheme utilizing 51 grid points within the domain — 4 < 77 < 4. 
The independency of the results to the selected number of grids and the size of the domain 
was confirmed. 


The procedure followed for model assessment in the spatially developing mixing layer is 
essentially the same, but the input moments are obtained by DNS. The computational 
procedure is based on a hybrid finite difference- pseudospectral scheme. A fifth-order compact 
parameter difference scheme (Carpenter, 1990) is used for discretization in the streamwise 
direction and a spectral collocation method (Givi and Madnia, 1993) employing Fourier 
expansions for the cross-stream discretization. Time discretization is achieved by the Adams- 
Bashforth scheme. The computational domain is bounded by 0 < x < 646 w ,— 166 w < 
y < 16<5u,, where 6 W is the vorticity thickness at the inflow. The resolution consists of 
256 finite difference points and 128 collocation nodes. With this resolution, reliable DNS 
with a Reynolds number Re = 200 and Damkohler numbers up to Da = 10 (based on the 
mean streamwise velocity difference and the vorticity thickness at the inlet) is possible. The 
inflow condition for the streamwise mean velocity profile is given by a hyperbolic tangent 
function. The free stream velocity ratio is set arbitrarily equal to = 3. The flow is forced 
randomly at the inflow to expedite the formation of large scale structures within the domain 
considered. The maximum amplitude of the perturbations is set equal to 5% of the mean 
flow. For a detail description of the numerical procedure we refer to Givi and Jou (1988); 
McMurtry and Givi (1992). 
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3.2 Frequency-Spectra of Reacting Scalars 


With the assumption of ergodic flow, the temporal correlations of the scalar mass fractions 
are defined (Bendat and Piersol, 1986). 


Raa(t) = Y A (t)Y A (t + r) (16) 

is the autocorrelation of the mass fraction of species A, and 

RAB{T) = Y A {t)Y B {t + T) (17) 

is the cross correlation of the mass fractions of the two species A and B. Associated with 
these correlations are the autocovariance function and the cross- variance function, defined 
respectively as: 


,-r 2 


C A a(t) = R A a(t) - Y a , Cab(t) = Rab(t) - Y A Y B . 


(18) 


The autospectral density function of scalar A is the Fourier transform of its autocorrelation: 


S AA (ft) = f Raa(t ) exp(—2xClTj)dT, j = I , (19) 

J — OO 

fl is the frequency. Similarly the cross-spectral density function for scalars A and B is: 


S A b(^) = / i? J 4 fl(r)exp(— 2jrflrj)dr. 

J — OO 

It is customary to express the cross-spectrum as: 


( 20 ) 
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Sab( ft) = 15^(0)! exp [-jQ^fi)] 


( 21 ) 


where 0(fi) is the “phase,” and the “coherence function” is defined by: 


\SabWI 

SaaWSbbW 

It is desirable to relate the time-correlations and the frequency spectral densities of the 
reacting scalars to the frequency-spectrum of the mixture fraction. This, in general, is a 
formidable task as the relation (if it could be obtained) is dependent on the chemistry. Since 
the primary influence of chemical reaction is at diffusion scales, it is impossible to provide 
simple mathematical relations amongst the various spectral density functions. In the limit 
of equilibrium chemistry, the situation simplifies considerably as shown by Kosaly (1993). 
However, an exact determination still requires multi- (in this case, two-) time level stochastic 
analyses. In the presence of spatial symmetry, it is possible to simplify the relations further. 
For example, at the center of the scalar mixing layer in the BSK experiments (tj = 0,f = |), 
all the statistical properties of the two reactants are identical. Based on this symmetry it is 
easy to show for Da —* oo (Kosaly, 1993): 




Raa{ t )U=0,D<i->oo = ■^BB( 7 ’)|i)=0,Da— »oo = “ [^Ff(t) + ^F*F*(' r )] TJ= o > (23) 

R-Ab{t ) |tj= 0 ,Da— >00 = RBA{r)\ n = 0 ,Da-*oo = ^ [~Rff(t) + Rf>F>(t)]„ =0 • (24) 

where F = 2/ — 1 and F" = \F\. Similar relations hold for the spectral densities of A 
and B in terms of those of F and F'. Therefore, for the spectral densities normalized 
by their respective variances, it is straightforward to show (from here on, the subscript 
r\ = 0,-Da — * oo is dropped): 
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£u(ft) = = -L* [s FF (si)F” + 5 F . F .(n)f*^ 

*A 4 /^ 


( 25 ) 


fraW 1 / S FF (ft)F* \ 

2 

To complete the formulation, the statistics of F, |F|, A and B are required. Kosaly (1993) 
assumes a gaussian pdf for F from which all the statistics are subsequently determined. This 
pdf is justified only at regions far away from the inlet. In this particular experiment this 
assumption is reasonable as the data show an approximate gaussian pdf (not considering the 
discrepancy associated with the infinite support). In other homogenous flows, however, this 
approximation can yield significant errors (Frankel et al., 1993) since at the initial stages of 
mixing the pdf is composed of two approximate delta functions, i.e. complete segregation. 
At such stages, the use of the PF generated pdfs is more reasonable. In particular, the 
Pearson Type II is justified at all the locations along the center of the layer. With this 
frequency, Eqs. (3)-(4) with a = /? yield: 


f2 — f>2 — f.2 — ^ 


2a + 1 


= £ 2 . 


(27) 


Therefore: 


= t - (!)’ • ^ - s ’ - G ’’ G < s > = Tmthj - m 

Substituting Eqs. (27)-(28) into Eqs. (25)-(26) yields: 

§aa(R) = S B b( fi) = 2^2 _ Q2 Sff(Q) + (29) 
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5 « (n ) = 1 L, (30) 

S AA (Q) I + &-CP 5f»f.(0) W 

S 5^jr(n) 

These equations would be identical to those suggested by Kosaly (1993) only when E — ► 0. 
In this limit, the Beta density approaches a gaussian pdf (Leemis, 1986); therefore: 


lim 

E— o 



7T — 2 


JT 


( 31 ) 


which is valid only at final stages of mixing in homogeneous flows. For the other limit- 
ing conditions near the inlet of the BSK’s mixing layer, or the initial stages of mixing in 
homogeneous flows (Madnia et al., 1992) we have E 2 = G = 1. Therefore, the physical 
requirement: 


&u(fi)Uo = (32) 

can be only realized via the relations provided by Eqs. (30)-(31). To demonstrate this point 
more clearly, the variations of the parameters appearing in Eqs. (29)-(30) are shown in Fig. 
2 in terms of E. In the form presented E = 1,0 denote the initial and the final stage of 
mixing, respectively. Note that only the values at E — ► 0 are in accord with Kosaly’s results. 

All the results presented here are based on the assumption of an asymptotic gaussian pdf for 
the variable F. However, the results of numerous recent laboratory measurements (Castaing 
et al., 1989; Sano et al., 1989; Castaing et al., 1990; Jayesh and Warhaft, 1991; Jayesh and 
Warhaft, 1992; Gollub et al., 1991; Lane et al., 1993) and numerical simulations (Metais and 
Lesieur, 1992; Kimuraand Kraichnan, 1993; Kerstein, 1991; Pumir et al., 1991) suggest that 
the exponential distribution may provide a more reasonable representation of the pdf in high 
Reynolds number flows. The family of exponential pdfs is represented by: 
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PF(ip,q) = C(K, q) exp (- - jp- - 

\ K{q,v) 


( 33 ) 


where a is the standard deviation and C = [2/f«r( 2 ^-)] -1 . The value q = 2 corresponds 
to a gaussian pdf and q = 1 implies a Laplace (double exponential) density. To use this 
density for the evaluation of the spectral density functions, Eqs. (25)-(26) are represented in 
the simpler form: 


*S/m(0) = Sbb{£1) = 


UTT + lTTi 


(34) 


£4b(A) _ J 2 

Saa(&) i -l 


(35) 


where U — F* n / F* 2 . Following the same procedure as outlined above, Eq. (33) yields: 


U = 1- 


r 2 (J) 

?r(f)r(*±i)' 


(36) 


This implies that the parameter U(q) is in the approximate range 0.25 < U < 1, depending 
on the degree of the exponential at the asymptotic stage of mixing. For q = 2, the value 
U = 0.36 is that suggested by Kosaly (1993). 

Based on the results presented here, we recommend Eqs. (29)-(30) for relating the autospec- 
tral and the cross-spectral density functions of the reacting scalars to the spectrum of the 
mixture fraction in homogeneous flows with symmetric pdf of the mixture fraction. These 
results are valid if the asymptotic stage of mixing yields gaussian statistics. For the general 
family of the exponential scalar pdfs, the corresponding relations are given by Eqs. (34)-(36). 
Again, it is emphasized that these relations (and also those provided by Kosaly (1993)) are 
only valid in the limit of equilibrium chemistry. It is impossible to generalize these results 
for non-equilibrium chemistry flows. Also, it is very difficult, perhaps impossible, to provide 
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analogous relations when the pdf is non-symmetric, i.e. at off-center regions of the labora- 
tory flow and/or in other homogeneous flows where the local compositional structure is not 
in stoichiometric proportions. 


4 Results 

4.1 Modeling of the Reactant Conversion 

The results of our parametric study indicates that the best overall match between the model 
predictions of the mixture fraction and the experimental data of BSK is established with 
Ci = 0.89 and C 2 = 5.72. LaRue et al. (1981) suggest C\ = 0.89, C2 = 2.25. The difference 
in the magnitude of C 2 is attributed to the observation that the peak of the scalar variance in 
the measurements of LaRue et al. (1981) is about 30% higher than that in BSK. In BSK, data 
are provided of reacting flows with several values of the Damkohler numbers, stoichiometric 
mixture fractions, and the Reynolds number. In all the cases, the measured results for the 
mean value of the mixture fraction agree well with model predictions as shown in Fig. 3(a). 
In this and subsequent figures 77* = tj/6, where 6 = rj(f = 0.9) — T)(f = 0.1). The measured 
variance of the mixture fraction show a dependence on the Reynolds number. The results in 
Fig. 3(b) indicate a reasonable agreement between the predicted results and the experimental 
data for Re = 5, 300. For Re = 11, 700, the model fails to capture the “valley” at the center 
of the layer. This is due to the gradient diffusion closure and also the independency of the 
model to the Reynolds number. Transverse variations of the skewness and kurtosis of the 
mixture fraction are shown in Figs. 4 and 5 for Re = 11,700. The model results are based 
on the Beta frequency yielding closed form analytical expressions for higher order moments. 
While the accordance of the predicted results with experimental data is encouraging, the 
local minima and maxima portrayed by the data are not well-predicted. Noticeably, the 
“kink” in the skewness profiles as suggested experimentally is not captured. The results 
based on this pdf for predicting the statistics of reacting scalars under equilibrium condition 
are given in Figs. 6 and 7. These results are obtained by closed-form analytical expressions 
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based on Eqs. (4)-(5). In Fig. 6, the cross-stream variations of the normalized mass fractions 
of the reactants are in reasonable accord with measured data. A somewhat similar level of 
agreement is observed in the profiles of the reactants’ standard deviations as Fig. 7 shows 
that both the distribution and the peak of the deviations are predicted well. The maximum 
deviation is near the center of the layer. 

Of greater interest is the statistical behavior of the reactants’ mass fractions in the non- 
equilibrium chemistry flow. On account of the fact that the statistics of mixture fraction 
should be invariant with the reaction rate, it is required to tailor the model accordingly. 
Therefore, in this case the magnitudes of the empirical constants C\ and C 2 are prescribed 
in such a way as to yield the best match to the first two measured moments of the reac- 
tants’ mass fractions, of course in accord with the respective parameterization scheme. The 
statistics of the mixture fraction are evaluated subsequently. This procedure is followed for 
both of the cases reported in the BSK experiments: Da = 0.3 and Da = 1.81. The com- 
parison between the predicted standard deviation of the mixture fraction and the measure 
data is shown in Fig. 8. With the degree of freedom available in the matching, the Dirichlet 
density parameterized by the scalar-energy yields better overall agreements. However, the 
results are not Damkohler number invariant. The consequence of this will be discussed later. 
Figures 9 and 10 show the profiles of skewness and kurtosis of the mixture fraction. Both 
Dirichlet closures perform well and the results are in a better agreement with data than 
those based on the univariate Beta density (Figs. 4-5). It is appropriate here to note that 
the multivariate gaussian pdf yields uniform values of higher order moments (skewness = 0, 
and kurtosis = 3); thus despite possessing more statistical information, i.e. a higher-degree 
of freedom for parameterization, this pdf is less superior to the Dirichlet density. The final 
performance test of the models is their capability of predicting the statistics of the reacting 
scalars under non-equilibrium conditions. In Figs. 11 and 12 results are presented of the 
mean product mass fraction and the scalar-energy, respectively. For consistency, the empir- 
ical constants are optimized for each model such that the predicted statistics of the mixture 
fraction are in accord with the experimental data. These figures also confirm the relative 
better performance of the scalar-energy parameterized Dirichlet. 


19 



A more elaborate validation of the models requires examination of the higher-order joint 
moments. Since these moments are not measured in the experiments, they are provided here 
by DNS of the spatially developing mixing layer. The structure of this layer as depicted 
by DNS is shown in Fig. 13, where a plot of the instantaneous mixture fraction contours is 
presented. Statistical analyses of the data are conducted at regions far away from the inflow 
so that the influence of random perturbations are sufficiently sensed. For statistical sampling, 
4,000 realizations are used. The sufficiency of this number of realizations was confirmed 
by detailed comparisons against results obtained at higher sampling rates. The conclusions 
drawn from the trends discussed below are independent of the streamwise location considered. 
Thus results are presented here only for x = 47.5£„, denoted by station I in Fig. 13. 

A detailed assessment of both Dirichlet closures indicates the superiority of the one parame- 
terized with the scalar-energy. For example, Figs. 14 and 15 show the cross-stream variation 
of the moment Y*Y B for two values of the Damkohler number. These results support the use 
of the scalar-energy for pdf parameterization. Also, the figures suggest that the agreement 
improves somewhat as the Damkohler number is increased. However, the use of neither of 
these models can be recommended when the chemistry is in complete equilibrium. In the 
limit of infinitely large Damkohler number all the joint moments of the scalar are identically 
zero. The use of the covariance Y'^Y' B — —YaYb yields = £ 2 = £3 = 0; thus the model fails. 
The use of the scalar-energy avoids this failure but is incapable of satisfying the orthogonal- 
ity condition (Eq. 4). In non-equilibrium flows, as the rate of reaction is increased but with 
finite Damkohler numbers, the covariance-quantified distribution quickly approaches the sin- 
gularity. This is demonstrated by a comparison of Figs. 14 and 15 where the results indicate 
a mild improvement in the scalar-energy parameterized model, but little or no change in the 
model based on the covariance. 

Despite the caveat on the use of Dirichlet density for flows under equilibrium chemistry, it 
is useful to examine its behavior in modeling of higher order moments. Figures 16-17 show 
the profiles of the joint central moments Y^Yg, and Y'^Yg. The results indicate that the 
Dirichlet multivariate density parameterized by the scalar-energy and the univariate Beta 
frequency yield a similar level of agreement with the DNS data. No distinct superiority is 
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established via the use of either of the closures. However, due to univariate nature of the 
Beta pdf, the orthogonality condition is identically satisfied. 

Finally, it is useful to compare the predicted statistics of the mixture fraction with the 
DNS data in non-equilibrium flows. This comparison is made in Figs. 18 and 19 for the 
skewness for both Damkohler numbers. The results are also compared with the predictions 
based on the univariate Beta pdf as parameterized by the DNS data. These figures indicate 
that both closures behave somewhat similarly without a dominant superiority of any of the 
two. This trend is observed for all the other moments of the mixture fraction (not shown 
here). However, the Dirichlet density does not yield an exact Beta frequency for the mixture 
fraction. It is straightforward to demonstrate that with the assumption of a Dirichlet density 
for the reacting scalars, the univariate pdf of the mixture fraction is a Beta density if the 
chemistry is frozen. That is, along the line V’i + = 1 for n — 2 in Eq. (7). The Dirichlet 

density yields a marginal Beta pdf for each of the scalars; thus a Beta distribution cannot 
be obtained for the mixture fraction. Therefore, in the limit of infinitely large Damkohler 
number, the use of the univariate Beta density is more appropriate. An optimum value of 
a finite Damkohler number below which the Dirichlet density is to be recommended, cannot 
be specified. Due to computational limitations, it is not possible to perform DNS of very 
large (but finite) Damkohler numbers. For the highest reaction rate considered in the BSK 
experiments, the Dirichlet density performs reasonably well. 


4.2 Frequency-Spectra of Reacting Scalars 

Many important spectral features of the reacting scalars can be exhibited by results extracted 
from the DNS of the spatially developing mixing layer. Here these results are presented of 
the layer under the conditions of frozen, equilibrium and non-equilibrium chemistry. The 
normalized autospectrum of scalar A is shown in Fig. 20 for several values of the Damkohler 
number at streamwise locations of x = 4 8 U and x = 166 w (stations II and III in Fig. 13) 
along the centerline ( y = 0). This figure indicates that at low frequencies, the influence of 
chemical reaction is not significant as the spectrum is influenced primarily by hydrodynamics. 
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At fl = 0 the amplitude must decrease as the Damkohler number is increased in order to 
accommodate for the decrease of the integral time-scale due to chemistry. An opposite trend 
is observed at the other end of the spectrum. At high frequencies, the autospectral density 
becomes more energetic as the magnitude of the Damkohler number is increased, and at 
the equilibrium limit the spectrum range extends to very high frequencies. This behavior is 
not observed in the results of BSK and Kosaly (1993) since a small frequency range, limited 
at the smallest scales to the hydrodynamics frequency is considered. In all the cases, the 
autospectrum exhibits local peaks at fi = 1.086Hz (for x = 4 8 W ) and at fi = 0.543Hz (for 
x = 16<5 W ). These values correspond approximately to the most unstable frequency of the 
mean flow and its first subharmonic (Michalke, 1965). Note that x = 4 S w corresponds to the 
approximate location where the first vortical roll-up occurs, and at x = 16^ the pairing of 
neighboring vortical structures is observed. After this point, the additional vortex pairings 
which occur downstream yield the shift of the spectra peaks towards higher frequencies. 
Since the area under the normalized autospectrum must be identical in all the cases (Bendat 
and Piersol, 1986), the increase of energies at high frequencies must be compensated by an 
energy decrease at low frequencies. This causes the curves of the spectral density functions 
to cross over each other at some intermediate values of the frequency. This behavior is 
observed in Fig. 20 and also in the results of BSK and Kosaly (1993). 

In Fig. 21 results are presented of the normalized autospectral density function of F , F*, and 
the reactant A under equilibrium condition. Although Saa and Sf>f- can be determined 
from the F-time series data, they cannot be calculated from Sff directly unless the exact 
two-time pdf of the F- process is known. Note that unlike the conditions at the center of 

A 

the symmetric configuration in the BSK experiments, here Sff • ^ 0. Thus, Eq. (29) or 
Eq. (34) cannot be used. However, the DNS results do indicate that the amplitude of Saa 
falls between Sff and Sf-f* at high frequencies. At low frequencies, due to the influence 
of hydrodynamics in causing asymmetric effects, the use of our analytical relations is not 
appropriate. These figures also show that the local peaks in the F- and A-autospectra axe 
not observed in the F*-autospectrum. This is expected since the spectrum of the absolute 
value of a process forced at a single particular frequency (like the process corresponding to 
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the transport of a single-eddy), contains energy at all the harmonics of that frequency. Since 
the absolute value of a process represents a non-linear system, the analytical determination 
of its spectrum is difficult. However, it is reasonable to expect that the autospectrum of 
F * decays slower than that of F at high frequencies and the results in Fig. 21 do show this 
behavior. 

The cross-spectral density function is characterized by the coherence function and the phase. 
High coherence between two random processes usually indicates that the two processes are 
either exactly in phase or out of phase (0 = 0, or 7r). Low coherence implies indeterminate 
phase values (Bendat and Piersol, 1986). In Figs. 22 and 23 results are presented of the 
coherence and the phase in the cross-spectral density functions of scalars A and B for different 
magnitudes of the Damkohler number. In the frozen chemistry limit, the coherence is unity 
and phase is equal to x at all frequencies. This is expected as the reactants are completely 
out of phase. At low frequencies the coherence is high and the non-equilibrium chemistry 
results are bounded by those corresponding to the frozen and equilibrium limits. In this 
range, the magnitudes of the phase for all cases are near x. At frequencies which correspond 
to the local peaks of the autospectra (Fig. 20), there are also high local coherence peaks. 
At the intermediate frequencies where the autospectra curves cross, the coherence drop and 
also cross over each other. The magnitude of the phase in this frequency range is random 
and as speculated by BSK this behavior is associated with the reactive-diffusive balance in 
the spectral transfer due the large-scale eddies and chemical reaction. After dropping to 
near zero at intermediate frequencies, the coherence rise at higher frequencies with a reverse 
effect of the Damkohler number. That is, the coherence increases as the Damkohler number 
increases. In this range, the phase is near zero (0 = 0, 2x) for all reacting cases. This 

A A A A 

behavior is in fact predicted by Eq. (35) in that as Sff 0, Sf’F* ^ 0 ( Sf»f • >> Sff), 
then Sab = Saa ■ Kosaly (1993) also recognizes the possibility of increasing coherence with 
a zero phase at high frequency. However, due to lack of experimental data this behavior 
could not be observed. The DNS results here capture this trend. 
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5 Concluding Remarks 


The experiments of Bilger et al. (1991) provide an extensive data set for assessing the 
role of turbulence fluctuations on the rate of reactant conversion and on the spectra of 
reactive scalars in turbulent shear flows. The objective here is to provide mathematical 
models to reproduce these data with the hope of suggesting simple working relations for 
predictive applications. Based on comparative assessment against these measured data and 
additional data generated here by direct numerical simulation, we recommend the Pearson 
family of probability density functions for statistical description of scalar transport in flows 
of this type. In particular, the Dirichlet frequency parameterized by the scalar-energy is 
recommended, in the absence of better alternatives, for modeling of the reactant conversion 
rate in non-premixed reacting mixing layers under non-equilibrium conditions. In the limit 
of frozen chemistry this density yields a Beta frequency for the marginal pdf of each of the 
reactants. This univariate frequency is recommended for stochastic treatment of both frozen 
and equilibrium flows. 

In the context considered, the models do not yield a consistent limiting condition for equilib- 
rium flows. In this limit, the Dirichlet density does not satisfy the orthogonality condition of 
the reactants and also does not yield a Beta density for the mixture fraction. This is due to 
inability of the distribution to include all the first and second order moments in its param- 
eterization. This problem is well-recognized in statistics and (classical) biometric literature 
(Johnson, 1949a). In principle, it is possible to construct a modified multivariate Pearson 
density which overcomes this predicament (Johnson, 1987). However, the parameters of the 
model cannot be algebraically related to input moments (that is, the pdf cannot be analyt- 
ically integrated); thus the model cannot be recommended for practical applications. The 
other known multivariate frequencies such as the joint gaussian (and distributions gener- 
ated by the JET and other schemes), do not share this problem but they do not possess 
appropriate physical properties to justify their use for combustion applications. 

The other alternative for pdf description is to utilize a transport equation governing its 
evolution (Pope, 1985). However, despite concentrated recent efforts it is not yet clear what 
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physical closure and/or computational schemes are to be used in the implementation of 
this scheme. Currently, the most promising methodology for this purpose is the Amplitude 
Mapping Closure (Chen et al ., 1989; Kraichnan, 1989), but it is not yet at a level suitable for 
modeling of multi-scalar transport in non- homogeneous flows (Pope, 1991). The results of 
extensive recent work (Miller et al., 1993; Frankel et al., 1993; Madnia et al ., 1992) indicate 
that at situations where the AMC can be enacted, other schemes based on the PF and the 
JET perform equally well. Amongst all available schemes, the PF is the simplest to use and 
until the shortcomings associated with the use of the pdf transport equation are resolved, 
the PF (or other parameterized pdfs prescribed based on known physics) will likely remain 
as the method of choice for practical applications. 

With the use of the PF family of pdfs, the spectral density functions of the reacting scalars are 
related to the frequency-spectrum of the mixture fraction. These results are more general 
than those provided previously in that the evolution of the mixture fraction pdf from an 
initial double-delta profile to an asymptotic gaussian distribution is taken into account. In 
the asymptotic limit of mixing completion, the results are in accord with those suggested 
by Kosaly (1993). The corrections to the frequency spectra for asymptotic exponential pdfs 
are also provided. These mathematical relations are valid for applications in homogeneous 
mixing layers such as the configuration considered by BSK. However, the trends portrayed 
by these relations are also in accord with those provided by DNS of the spatially developing 
reacting mixing layer. 

At this point it is instructive to reiterate the limitations of the model in the format utilized 
in this work. The primary drawbacks are associated with: (1) modeling of turbulent fluxes, 
(2) the single-point nature of the statistical description (3) the restrictive applicability of 
the equations governing the spectral densities, and (4) the non-generality of the frequencies. 
The first three problems are also encountered in approaches based on a pdf transport; the 
last problem is particular to the model here. The first problem manifests itself in the lim- 
itations of the gradient diffusion closure and is present in any pdf approach in which the 
statistics of the velocity field are not included (Givi and McMurtry, 1988). It is now well- 
recognized that this closure is not capable of capturing the spatial transport of the scalar 
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pdf in non-homogeneous turbulence, particularly when the flow is dominated by organized 
coherent structures (Koochesfahani and Dimotakis, 1986). To remedy this problem, one can 
potentially use higher order models such as spectral closures (Frankel et al., 1992); but the 
implementation of these models in non-homogeneous flows is extremely difficult, if not im- 
possible. The second problem infers that the errors associated with modeling of the first two 
moments of the scalar result in the contamination of the pdf and thus all the higher order mo- 
ments. To rectify the situation, one can use multi-point statistical closures; but again their 
use is not practical at this point (Pope, 1990). The third limitation is understandable as the 
spectral manipulation of nonlinear systems is a challenging task, which sometimes is not even 
recommended (Bendat, 1990). Here with the assumption of equilibrium chemistry and with 
the spatial symmetry of the problem, it is possible to relate the spectra of the reactants to 
the frequency spectrum of the mixture fraction. For problems without such a symmetry, the 
final results can be only obtained by the numerical integration of the equations relating the 
correlation functions. In non-equilibrium flows it is impossible to provide analogous relations 
due to closure problems encountered in determining the correlation functions (Corrsin, 1958; 
Lee, 1966; Corrsin, 1981). Even in the form presented, the final equations require the knowl- 
edge of the spectral densities of the absolute value of the mixture fraction. Without pertinent 
available data (from laboratory experiments or DNS), the evaluation of these densities require 
the knowledge of “two-time process” of the signal. Note that for the simple case of an ergodic 
gaussian stochastic process, an analytical evolution of these spectral densities is infeasible 
as realized in the literature on nonlinear system identifications (Bendat and Piersol, 1986; 
Bendat, 1990). Finally, the drawbacks associated with the non-generality of the present 
model, or other approaches based on parameterized pdfs, are well-recognized. In applying 
the PF generated pdf for modeling of other types of turbulent combustion systems ( e.g pre- 
mixed flames, temperature dependent reaction rates, etc.) one must have some knowledge 
of the systems before utilizing a priori schemes. 
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Figure captions 


Figure 1. Schematic diagrams of (a) the homogenous scalar mixing layer, and (b) the spatially 
developing mixing layer. 

Figure 2. Variation of the parameters of Eqs. (29)-(30). The points identified by • correspond 
to the results of Kosaly (1993). 

Figure 3(a). Cross-stream variation of the mean value of the mixture. Frozen chemistry. 

Figure 3(b). Cross-stream variation of the standard deviation of the mixture fraction. Frozen 
chemistry. 

Figure 4. Cross-stream variation of the skewness of the mixture fraction. Frozen chemistry. 

Figure 5. Cross-stream variation of the kurtosis of the mixture fraction. Frozen chemistry. 

Figure 6. Cross-stream variation of the mean values of the reactants’ mass fractions under 
equilibrium chemistry. 

Figure 7. Cross-stream variation of the standard deviations of the reactants’ mass fractions 
under equilibrium chemistry. 

Figure 8. Cross-stream variation of the standard deviation of the mixture fraction. 

Figure 9. Cross-stream variation of the skewness of the mixture fraction. 

Figure 10. Cross-stream variation of the kurtosis of the mixture fraction. 

Figure 11. Cross-stream variation of the mean product mass fraction. 

Figure 12. Cross-stream variation of the scalar-energy. 

Figure 13. Plot of the mixture fraction contours. 

Figure 14. Cross-stream variation of the moment Y^Yg 2 for Da = 0.3. 

Figure 15. Cross-stream variation of the moment Y^Yg 2 for Da = 10. 

Figure 16. Cross-stream variation of the moment Yj^Yg 2 . Equilibrium chemistry. 

Figure 17. Cross-stream variation of the moment YfYb. Equilibrium chemistry. 

Figure 18. Cross-stream variation of the skewness of the mixture fraction (Da = 0.3). 

Figure 19. Cross-stream variation of the skewness of the mixture fraction (Da = 10). 

Figure 20. Autospectrum density function of the mass fraction of species A at (a) x = 
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4<5 w ,y = 0, (b) x = 16<5 W , y = 0. 

Figure 21. Autospectra of F, F m and A (Da — ► oo) at (a) x = y = 0, (b) x = 165 w , y = 0. 
Figure 22. Coherence between reactants A and B at (a) z = 4 8 w ,y = 0, (b) z = 166 w , y = 0. 
Figure 23. Phase between reactants A and B at (a) z = 4 8 w ,y = 0, (b) z = 16^, y = 0. 
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